set more off
clear all

version 17.0



global ylist d_wood d_charcoal d_lpg lq_coal_hat d_qcoal lq_coal_hat2 lexp_lpg d_elpg lexp_lpg2
global ylist2 d_wood d_charcoal d_lpg lq_coal_hat d_qcoal lq_coal_hat2 totexp_lpg_year d_elpg totexp_lpg_year2
global xlist linccap lsize housing_nrooms housing_selfowned sexhead lagehead lyeduc agrihead unemphead cook_jobweeks cook_female lyeduccook elec d_road
global xlist2 linccap lsize housing_nrooms housing_selfowned sexhead lagehead lyeduc agrihead unemphead cook_jobweeks cook_female lyeduccook rural elec d_road

************************************************************************************************
*** 7	Welfare, health and environmental impacts ***************************************************************
************************************************************************************************

*** mlogit for eliciting differences in firewood and charcoal consumption between hh's mainly using firewood and those mainly using charcoal
use "$datain\glss_combined.dta", clear
keep if glss==6 
keep if cooking<=3 & cooking>=1 
egen district_group=group(region district) 

replace d_road=1 if urban==1
replace elec=1 if urban==1
replace d_road=. if d_road==99
replace elec=. if elec==99

tab ethnicgroup, g(deth)
tab religion, g(drel)

set matsize 2000
/*activate this code to see kg-estimates 
mlogit cooking $xlist2 drel2 drel3 drel4 drel5 drel6 drel7 drel8 deth2 deth3 deth4 deth5 deth6 deth7 deth8 deth9 i.district_group , baseoutcome(1)
predict p
gen tag=(p>=0.9999|p<=0.0001)
drop p		

mlogit cooking  $xlist2 drel2 drel3 drel4 drel5 drel6 drel7 drel8 deth2 deth3 deth4 deth5 deth6 deth7 deth8 deth9 i.district_group  if tag==0 , baseoutcome(1) 
est sto mlog
drop tag

predict pwood if e(sample)
predict yhat_1, outcome(1)
predict yhat_2, outcome(2)
predict yhat_3, outcome(3)

gen qcoal_year=q_coal_hat*12
bysort urban: sum qcoal_year if treat==0 & cooking==1 & yhat_2>=0.5 //consumption of main firewood users with high probability of mainly using charcoal
bysort urban: sum qcoal_year if treat==0 & cooking==2 & yhat_1>=0.5 //consumption of main charcoal users with high probability of mainly using firewood
*/

*** further statstics for computations further below ***
use "$datain\glss_combined.dta", clear
keep if glss==6 
bysort urban: sum d_wood d_charcoal d_lpg if treat==0 & cooking!=0 & cooking!=.
bysort urban: sum hhsize n_child_u5
sum wta_s
egen temp=total(wta_s)
drop temp
bysort urban: egen temp=total(wta_s)
tab temp
drop temp

tab cooking if treat==0 & cooking!=0
gen totexp_charcoal_year=totexp_charcoal*12
gen totexp_wood_year=totexp_wood*12
replace totexp_Diesel=totexp_Diesel*12
replace totexp_Petrol=totexp_Petrol*12

foreach var in totexp_wood_year totexp_charcoal_year totexp_lpg_year totexp_Diesel totexp_Petrol{
	replace `var'=`var'*(1/1.954) //to USD
	winsor2 `var', trim replace c(1 99)
}

bysort urban cooking: sum totexp_wood_year totexp_charcoal_year totexp_lpg_year if treat==0 & cooking!=0
bysort urban: sum totexp_wood_year totexp_charcoal_year totexp_lpg_year if treat==0 & cooking!=0

bysort urban: sum totexp_Diesel totexp_Petrol if treat==0 & cooking!=0


gen temp=totexp_wood_year>0 & totexp_wood_year!=.
bysort urban: sum temp if treat==0 & cooking!=0
sum temp if treat==0 & cooking!=0 & d_wood==1


gen cooktime_hr_month=cooktime_min*4/60
gen fetchwood_hr_month=fetchwood_min*4/60

bysort urban: sum eqsc hhsize n_child_u5 d_wood d_charcoal d_lpg 

gen qcoal_year=q_coal_hat*12

bysort urban cooking: sum qcoal_year if treat==0 & cooking!=0
bysort urban: sum qcoal_year if treat==0 & cooking!=0

*** time costs ***
gen time_wood_month=cooktime_hr_month+fetchwood_hr_month
bysort urban: sum cooktime_hr_month fetchwood_hr_month time_wood if d_wood==1
bysort urban: sum cooktime_hr_month if d_charcoal==1
bysort urban: sum cooktime_hr_month if d_lpg==1

collapse (sum) wta_s_hhsize wta_s (mean) totexp_wood_year totexp_charcoal_year totexp_lpg_year, by(urban)
gen perc_lower=(0.033-(1.96*0.011))*wta_s if urban==1
gen perc_upper=(0.033+(1.96*0.011))*wta_s if urban==1
replace perc_lower=0*wta_s if urban==0
replace perc_upper=(0.052+(1.96*0.027))*wta_s if urban==0

gen minwage=	5.24 			*(1/1.954)
*				minimum wage 	exchange rate
gen time_wood0=2*minwage*perc_lower*12/1000000 if urban==1 //annual time costs due to switching to fuelwood, assuming two days per month additional cooking time (difference to mean time spent when using charcoal)
gen time_wood1=2*minwage*perc_upper*12/1000000 if urban==1 
gen time_coal0=	(0.159-1.96*0.052)		*0.5 * minwage * 12 * wta_s/1000000 if urban==1 //annual time costs due to switching to fuelwood, assuming half a day per month additional cooking time (difference to mean time spent when using LPG)
gen time_coal1=	(0.159+1.96*0.052)		*0.5 * minwage * 12 * wta_s/1000000 if urban==1 


replace time_wood0=3*minwage*perc_lower*12/1000000 if urban==0 
replace time_wood1=3*minwage*perc_upper*12/1000000 if urban==0 

egen time_lower=rowtotal(time_wood0 time_coal0)
egen time_upper=rowtotal(time_wood1 time_coal1)
sum time_lower time_upper

egen time_total_lower=total(time_lower)
egen time_total_upper=total(time_upper)

********************************************************************************
sum time_to* //time costs, lower and upper bound
********************************************************************************

*** cooking fuel costs ***
gen fcost_wood0=(3.5					*1.2			*wta_s 				+  perc_lower	*(9.2-3)*1.2)/1000000	if urban==1
*				avg. wood fuel costs	price incr. 	incr. in # of hhs					diff. in wood exp. between hh's mainly using firewood & mainly using charcoal
gen fcost_wood1=(3.5*1.2*wta_s +  perc_upper*(9.2-3)*1.2)/1000000	if urban==1

replace fcost_wood1=(perc_lower*(3.19-3.05))/1000000	if urban==0 //assume negligible transport price increase for rural areas
replace fcost_wood0=(perc_upper*(3.19-3.05))/1000000	if urban==0 //negative, hence lower bound

gen fcost_coal0=(52.1*1.081*1.2-61.4) *wta_s/1000000 if urban==1
gen fcost_coal1=(52.1*1.275*1.2-61.4) *wta_s/1000000 if urban==1

replace fcost_coal1=0 if urban==0
replace fcost_coal0=(19.9*(1-.26616)*1.2-19.9) *wta_s/1000000 if urban==0 //negative, hence lower bound

gen fcost_lpg1=(15*.3)*wta_s/1000000	if urban==1 
replace fcost_lpg1=(2.3*.3)*wta_s/1000000 if urban==0 
gen fcost_lpg0=0 


gen fcost_tot0=(fcost_coal0+fcost_wood0+fcost_lpg0)
gen fcost_tot1=(fcost_coal1+fcost_wood1+fcost_lpg1)

egen fcost_total_lower=total(fcost_tot0)
egen fcost_total_upper=total(fcost_tot1)


********************************************************************************
sum fcost_total* //cooking fuel costs, lower and upper bound
********************************************************************************

*** transport fuel costs ***
gen fcost_transp_diesel=0.2*56.4*wta_s if urban==1
gen fcost_transp_petrol=0.2*0*wta_s if urban==1

replace fcost_transp_diesel=0.2*35.6*wta_s if urban==0
replace fcost_transp_petrol=0.2*0*wta_s if urban==0

gen fcost_transp=fcost_transp_diesel+fcost_transp_petrol
egen fcost_transp_total=total(fcost_transp)
replace fcost_transp_total=fcost_transp_total/1000000

********************************************************************************
sum fcost_transp_total //transport fuel costs
********************************************************************************


*** ghg emissions ***
//ghg emissions by fuel, emission factors from Bailis et al. 2015, three stone for wood, metal jiko for coal use, and Kenyan earth-mound kiln for charcoal prduction
gen ghg_wood_kg_sus =	(1584 * 1 *0.11		///CO2, adjusted by sustainable fraction of woodfuel harvest
					+57.4 * 3		///CO 
					+2.0 * 28		///CH4
					+7.8 * 12		///NMHC  
					+0.9 * 650		///BC 
					+1.5 * (-75))	///OC 
					/1000 			//in kg
label var ghg_wood_kg_sus "est. ghg emissions of 1 kg of woodfuels (kgCO2e)"

gen ghg_coal_use_kg_sus = 	(2857 * 1	*0.11	///CO2, adjusted by sustainable fraction of woodfuel harvest
						+195.9 * 3		///CO
						+8.9 * 28		///CH4
						+20.3 * 12		///NMHC
						+0.2 * 650		///BC
						+1.5 * (-75))	///OC
						/1000			//in kg
gen ghg_coal_pyrolysis_kg_sus = (1802 * 1	*0.11 ///CO2, adjusted by sustainable fraction of woodfuel harvest
							+223 * 3		///CO
							+44.6 * 28		///CH4
							+93 * 12		///NMHC
							+0.1 * 0		///NOx assumed!!
							+0.2 * 265)		///N2O
							/1000			//in kg
gen ghg_coal_kg_sus = ghg_coal_use_kg_sus+ghg_coal_pyrolysis_kg_sus
label var ghg_coal_kg_sus "est. ghg emissions of 1 kg of charcoal (kgCO2e)"

*what Net calorific value of LPG? (TJ/kg) https://www.seai.ie/data-and-insights/seai-statistics/conversion-factors/
gen NCV_lpg=0.047156/1000

*kg of co2e per TJ https://www.ipcc-nggip.iges.or.jp/public/2006gl/pdf/2_Volume2/V2_2_Ch2_Stationary_Combustion.pdf
gen ghg_lpg_kg=	63100*NCV_lpg	+5*28*NCV_lpg	+0.1*265*NCV_lpg
gen ghg_lpg_kg_sus=	63100*NCV_lpg	+5*28*NCV_lpg	+0.1*265*NCV_lpg
*				CO2				CH4				N2o
label var ghg_lpg_kg "est. ghg emissions of 1 kg of LPG (kgCO2e)"


gen MJ_coal=29.8 //The net calorific value estimates used for conversion of mass to energy are the values used in the Water Boiling Test 4.2.2 (WBT Technical Committee, 2013), which for charcoal is 29.8 MJ/kg, for wood is 19.0 MJ/kg, for LPG is 44.7 MJ/kg, and for kerosene is 39.7 MJ/kg.
gen MJ_wood=19
gen MJ_lpg=44.7


gen wood_diff_lower=(perc_lower*((	281.5-112		)*MJ_coal/MJ_wood)) if urban==1
*									diff in annual charcoal consumption between main firewood and main charcoal users
*									262.3-122.6 <-- avg. values
gen wood_diff_upper=(perc_upper*((281.5-112)*MJ_coal/MJ_wood)) if urban==1 
replace wood_diff_lower=0 if urban==0 
*replace wood_diff_upper=(perc_upper*((201.7-52)*MJ_coal/MJ_wood)) if urban==0 
*									  207.3-30.6 <-- avg. values


gen coal_diff_lower=(	177.42)									*(0.159-1.96*0.052)*wta_s if urban==1 
*						initial yrly avg. charcoal consumption	lower bound impact estiamte
gen coal_diff_upper=(177.42)*(0.159+1.96*0.052)*wta_s if urban==1 
replace coal_diff_lower=0 if urban==0  
replace coal_diff_upper=(((63.62)*(1-0.26616))-63.62)*wta_s if urban==0 
replace wood_diff_upper=-(coal_diff_upper*MJ_coal)/MJ_wood if urban==0 


gen lpg_diff_lower=-(coal_diff_lower*MJ_coal+wood_diff_lower*MJ_wood)/MJ_lpg if urban==1 
gen lpg_diff_upper=-(coal_diff_upper*MJ_coal+wood_diff_upper*MJ_wood)/MJ_lpg if urban==1
replace lpg_diff_lower=0 if urban==0 
*replace lpg_diff_upper=-(coal_diff_upper*MJ_coal+wood_diff_upper*MJ_wood)/MJ_lpg if urban==0
replace lpg_diff_upper=0 if urban==0

gen temp=lpg_diff_upper*MJ_lpg+wood_diff_upper*MJ_wood+coal_diff_upper*MJ_coal
sum wood_diff_* coal_diff_* lpg_diff_* temp if urban==0
sum wood_diff_* coal_diff_* lpg_diff_* temp if urban==1



foreach var in wood coal lpg{
	foreach case in lower upper{
		gen ghgdiff_`var'_sus_`case'=`var'_diff_`case'*ghg_`var'_kg_sus/1000/1000000 //in million tonnes
	}
}
sum ghgdiff*

egen yearly_lower_sus=rowtotal(ghgdiff_wood_sus_lower ghgdiff_coal_sus_lower ghgdiff_lpg_sus_lower)
egen yearly_upper_sus=rowtotal(ghgdiff_wood_sus_upper ghgdiff_coal_sus_upper ghgdiff_lpg_sus_upper)

egen ghg_lower_tot=total(yearly_lower_sus)
egen ghg_upper_tot=total(yearly_upper_sus)

gen annual_emissions=35340*1000/1000000 //emissions in 2012 (https://data.worldbank.org/indicator/EN.ATM.GHGT.KT.CE?locations=GH)

gen share_sus_upper=yearly_upper_sus/annual_emissions
gen share_sus_lower=yearly_lower_sus/annual_emissions

gen kg_woodextraction_lower=wood_diff_lower+coal_diff_lower*6 //assuming 6kg firewood per 1kg of charcoal
gen kg_woodextraction_upper=wood_diff_upper+coal_diff_upper*6
*sum kg_woode*

gen woodextraction_yrly_Mton_lower=(kg_woodextraction_lower/1000)/1000000
gen woodextraction_yrly_Mton_upper=(kg_woodextraction_upper/1000) /1000000

egen wood_extraction_lower_tot=total(woodextraction_yrly_Mton_lower)
egen wood_extraction_upper_tot=total(woodextraction_yrly_Mton_upper)

********************************************************************************
//emission totals
*sum share*  //emissions as share
*sum ghg_lower_tot ghg_upper_tot //total emissions
*sum wood_e* //additional wood demand

//by rural/urban
bysort urban: sum share* 
bysort urban: sum yearly_upper_sus yearly_lower_sus
bysort urban: sum woode*
********************************************************************************

*** health costs ***
import delimited "$datain\IHME_CRCurve_parameters.csv", clear 

gen MR=15.4 if cause=="ALRI" //mortality rates from Jeuland et al. (2018)
replace MR=1.81 if cause=="COPD"
replace MR=10.7 if cause=="IHD"
replace MR=1.27 if cause=="LC"

gen IR=0.638 if cause=="ALRI" //incidence rates from Jeuland et al. (2018)
replace IR=4.1 if cause=="COPD"
replace IR=1 if cause=="IHD"
replace IR=0.014 if cause=="LC"

replace MR=MR/10000
replace IR=IR/10000

gen size=3.652653 //average urban hh-size, from GLSS6
replace size=.4330423  if cause=="ALRI" //average number of children under 5 in ubran hh's, from GLSS6
gen size_r=4.752761 //rural
replace size_r=.6817841 if cause=="ALRI" //rural


keep if cause=="ALRI" | cause=="COPD" | cause=="IHD" | cause=="LC"
keep if age=="AllAge" | age=="All Age"

gen RR_wood=1+alpha*(1-exp(-beta*(423.95-zcf)^delta)) //firewood 
gen RR_coal=1+alpha*(1-exp(-beta*(156.16-zcf)^delta)) //charcoal
gen RR_lpg=1+alpha*(1-exp(-beta*(42.7-zcf)^delta)) //LPG

gen PAFu_case1=(((.169+0.11)*(RR_wood-1)+(.488-0.011+0.01)*(RR_coal-1)+(.343-0.01)*(RR_lpg-1))/(1+(.169+0.11)*(RR_wood-1)+(.488-0.011+0.01)*(RR_coal-1)+(.343-0.01)*(RR_lpg-1)))-(.169*(RR_wood-1)+.488*(RR_coal-1)+.343*(RR_lpg-1))/(1+.169*(RR_wood-1)+.488*(RR_coal-1)+.343*(RR_lpg-1))
gen PAFu_case2=(((.169+0.05456)*(RR_wood-1)+(.488-0.05456+0.03)*(RR_coal-1)+(.343-0.03)*(RR_lpg-1))/(1+(.169+0.05456)*(RR_wood-1)+(.488-0.05456+0.03)*(RR_coal-1)+(.343-0.03)*(RR_lpg-1)))-(.169*(RR_wood-1)+.488*(RR_coal-1)+.343*(RR_lpg-1))/(1+.169*(RR_wood-1)+.488*(RR_coal-1)+.343*(RR_lpg-1))

gen PAFr_case1=(((.781+0.10492)*(RR_wood-1)+(.166-0.09392)*(RR_coal-1)+(.053-0.011)*(RR_lpg-1))/(1+(.781+0.10492)*(RR_wood-1)+(.166-0.09392)*(RR_coal-1)+.053*(RR_lpg-1)))-((.781*(RR_wood-1)+.166*(RR_coal-1)+.053*(RR_lpg-1))/(1+.781*(RR_wood-1)+.166*(RR_coal-1)+.053*(RR_lpg-1)))
gen PAFr_case2=((.781*(RR_wood-1)+.166*(RR_coal-1)+.053*(RR_lpg-1))/(1+.781*(RR_wood-1)+.166*(RR_coal-1)+.053*(RR_lpg-1)))-((.781*(RR_wood-1)+.166*(RR_coal-1)+.053*(RR_lpg-1))/(1+.781*(RR_wood-1)+.166*(RR_coal-1)+.053*(RR_lpg-1)))

bysort cause: sum PAF*

gen Mortu1=size*(PAFu_case1)*MR
gen Mortu2=size*(PAFu_case2)*MR
gen Mortr1=size*(PAFr_case1)*MR
gen Mortr2=size*(PAFr_case2)*MR

gen Morbu1=size*(PAFu_case1)*IR
gen Morbu2=size*(PAFu_case2)*IR
gen Morbr1=size*(PAFr_case1)*IR
gen Morbr2=size*(PAFr_case2)*IR

egen Mortu_lower=rowmin(Mortu1 Mortu2)
egen Mortu_upper=rowmax(Mortu1 Mortu2)
egen Mortr_lower=rowmin(Mortr1 Mortr2)
egen Mortr_upper=rowmax(Mortr1 Mortr2)

egen Morbu_lower=rowmin(Morbu1 Morbu2)
egen Morbu_upper=rowmax(Morbu1 Morbu2)
egen Morbr_lower=rowmin(Morbr1 Morbr2)
egen Morbr_upper=rowmax(Morbr1 Morbr2)
sum Mor*

gen VSL=30000 //value of statistical life from Jeuland et al. (2018)

gen COI=33 if cause=="ALRI"
replace COI=67.5 if cause=="COPD"
replace COI=45 if cause=="IHD"
replace COI=1650 if cause=="LC"

gen cl1=0.7 if cause=="ALRI"
replace cl1=0.3 if cause=="COPD"
replace cl1=0.2 if cause=="IHD"
replace cl1=0.2 if cause=="LC"

gen cl2=0.1 if cause=="ALRI"
replace cl2=0.2 if cause=="COPD"
replace cl2=0.1 if cause=="IHD"
replace cl2=0.1 if cause=="LC"

gen cl3=0.2/3 if cause=="ALRI"
replace cl3=0.5/3 if cause=="COPD"
replace cl3=0.7/3 if cause=="IHD"
replace cl3=0.7/3 if cause=="LC"

gen cl4=0.2/3 if cause=="ALRI"
replace cl4=0.5/3 if cause=="COPD"
replace cl4=0.7/3 if cause=="IHD"
replace cl4=0.7/3 if cause=="LC"

gen cl5=0.2/3 if cause=="ALRI"
replace cl5=0.5/3 if cause=="COPD"
replace cl5=0.7/3 if cause=="IHD"
replace cl5=0.7/3 if cause=="LC"

//without discount rate
/*gen Value_Mort_lower=(Mort_lower)*VSL/12
gen Value_Mort_upper=(1+pi)*(Mort_upper)*VSL/12
gen Value_Morb_lower=(1+pi)*(Morb_lower)*COI/12
gen Value_Morb_upper=(1+pi)*(Morb_upper)*COI/12*/

//with discount rate
gen Value_Mortu_lower2=((Mortu_lower*cl1/(1))+(Mortu_lower*cl2/(1+delta))+(Mortu_lower*cl3/(1+delta)^2)+(Mortu_lower*cl4/(1+delta)^3)+(Mortu_lower*cl5/(1+delta)^4))*VSL/12
gen Value_Mortu_upper2=((Mortu_upper*cl1/(1))+(Mortu_upper*cl2/(1+delta))+(Mortu_upper*cl3/(1+delta)^2)+(Mortu_upper*cl4/(1+delta)^3)+(Mortu_upper*cl5/(1+delta)^4))*VSL/12
gen Value_Morbu_lower2=((Morbu_lower*cl1/(1))+(Morbu_lower*cl2/(1+delta))+(Morbu_lower*cl3/(1+delta)^2)+(Morbu_lower*cl4/(1+delta)^3)+(Morbu_lower*cl5/(1+delta)^4))*VSL/12
gen Value_Morbu_upper2=((Morbu_upper*cl1/(1))+(Morbu_upper*cl2/(1+delta))+(Morbu_upper*cl3/(1+delta)^2)+(Morbu_upper*cl4/(1+delta)^3)+(Morbu_upper*cl5/(1+delta)^4))*VSL/12

gen Value_Mortr_lower2=((Mortr_lower*cl1/(1))+(Mortr_lower*cl2/(1+delta))+(Mortr_lower*cl3/(1+delta)^2)+(Mortr_lower*cl4/(1+delta)^3)+(Mortr_lower*cl5/(1+delta)^4))*VSL/12
gen Value_Mortr_upper2=((Mortr_upper*cl1/(1))+(Mortr_upper*cl2/(1+delta))+(Mortr_upper*cl3/(1+delta)^2)+(Mortr_upper*cl4/(1+delta)^3)+(Mortr_upper*cl5/(1+delta)^4))*VSL/12
gen Value_Morbr_lower2=((Morbr_lower*cl1/(1))+(Morbr_lower*cl2/(1+delta))+(Morbr_lower*cl3/(1+delta)^2)+(Morbr_lower*cl4/(1+delta)^3)+(Morbr_lower*cl5/(1+delta)^4))*VSL/12
gen Value_Morbr_upper2=((Morbr_upper*cl1/(1))+(Morbr_upper*cl2/(1+delta))+(Morbr_upper*cl3/(1+delta)^2)+(Morbr_upper*cl4/(1+delta)^3)+(Morbr_upper*cl5/(1+delta)^4))*VSL/12

bysort cause: sum V*

collapse (mean) Value* , by(cause)
collapse (sum) Value*
gen pi=0.112
gen nhh=3656461
gen nhhr=2945023

gen Mortu_lower=(1+pi)*Value_Mortu_lower2 //ambient air pollution adjustment
gen Mortu_upper=(1+pi)*Value_Mortu_upper2
gen Morbu_lower=(1+pi)*Value_Morbu_lower2
gen Morbu_upper=(1+pi)*Value_Morbu_upper2
gen Mortr_lower=(1+pi)*Value_Mortr_lower2
gen Mortr_upper=(1+pi)*Value_Mortr_upper2
gen Morbr_lower=(1+pi)*Value_Morbr_lower2
gen Morbr_upper=(1+pi)*Value_Morbr_upper2

gen healthu_upper=(Mortu_upper+Morbu_upper)*nhh/1000000 //# urban hh's: 3656461
gen healthu_lower=(Mortu_lower+Morbu_lower)*nhh/1000000 //# rural hh's: 2945023
gen healthr_upper=(Mortr_upper+Morbr_upper)*nhhr/1000000 //# urban hh's: 3656461
gen healthr_lower=(Mortr_lower+Morbr_lower)*nhhr/1000000 //# rural hh's: 2945023

gen health_upper=healthu_upper+healthr_upper
gen health_lower=healthu_lower+healthr_lower

********************************************************************************
sum health_* //health costs
********************************************************************************


